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Coherent multiple scatterings of cc quark pairs in the environment of heavy ion collisions have been 
used in a previous work by Qiu et al. 0| to study J/tp suppression. That model suggests that heavy 
quark re-scatterings in a cold nuclear medium can completely explain the centrality dependence of 
the observed J/tp suppression in Pb+Pb collisions at the SPS ;2]. Their calculations also revealed 
significant differences under the assurnptions of a color singlet or color octet production mechanism. 
A more recent analytic calculation Q, which includes incoherent final-state re-scatterings with 
explicit momentum transfer fluctuations in three dimensions, indicates much less suppression and 
little sensitivity to the production mechanism. In this article, we study simultaneously both the 
J/ip suppression and pr modifications, at SPS and RHIC energies. We mainly focus on incoherent 
momentum transfer fluctuations in two dimensions, which is more appropriate for the heavy-ion 
collision kinematics. Our analytic and Monte-Carlo calculations reinforce the analytic results in 
Additionally, we find that the experimental J/ip suppression and (py) from nucleus-nucleus collisions 
at the SPS or RHIC cannot simultaneously be described in this incoherent multiple scattering 
framework for any value of the fluctuation strength parameter {k^). 

PACS numbers: 25.75.-q, 25.75.Nq 

I. INTRODUCTION 

Heavy quarkonia (cc, hh) have long been of great theoretical and experimental interest as sensitive probes of color 
deconfinement and thermalization in heavy-ion collisions 0,0]. The survival probability of these bound states depends 
on the density and (effective) temperature of the system, leading to the expectation of progressively suppressed 
quarkonia yields with increasing collision centrality and/or center-of-mass energy. 

The most extensively studied quarkonium state is the J/ip. The pedagogical picture of J formation in nucleus- 
nucleus collisions proceeds in multiple stages. First, two Lorentz contracted nuclei pass though one another and a 
particular partonic hard scattering forms a cc pair, a process requiring tec ^ ^ ^-^^ fm/c. The pair is then 
swept through the remaining fast traveling cold nuclear material, of length L, as shown in Figure ^ This crossing 
time is up to tcross ^ Diameternucieus / "t ~ 0.1 fm/c at RHIC (in the center-of-mass frame). A surviving cc pair can 
form a J/tp at ij/^ ^ Radius j/^/c ~ 0.3 fm/c unless the hot dense medium left in the wake of the nuclear collision, 
lasting tmedium ^ 10 fm/c, interferes. With a suitable model, J/^ measurements can therefore help extract, or at 
least constrain, the properties of the dense medium created in heavy-ion collisions. 

A theoretical framework for calculating J /ip suppression has been proposed by Qiu, Vary, and Zhang In their 
approach (QVZ), J/ip production factorizes into two stages: first perturbative production of cc pairs, followed at a 
much later stage by J /ip formation. The corresponding A + B ^ J /ip + X production cross section at leading-order 
in the strong coupling is dependent on the transition probability, F{q'^) for cc to evolve into a final J/4', where is 
the square of the relative momentum of the cc pair. Various J/ip formation mechanisms can be accommodated via 
different F functions. In QVZ, the cc J /ip transition probabilities in the color singlet and color octet channels were 
parameterized as 

= A^i/i^('7')exp[-gV(24)] (color singlet) (1) 

F'^lj/^{q') = N^^l6{q')e{Am'^-Arnl~q^)^[l-^-j^^y (color octet) (2) 

where TOc is the charm quark mass, m' is the mass scale for the open charm threshold {D meson mass), and Nj/^ 
and a f are parameters fixed from hadron-hadron collision data. 
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FIG. 1: Diagram of a nucleus-nucleus collision. Arrows at the center of each nucleus indicate the direction of travel. The dot 
in the center represents a nucleon-nucleon collision, and the distance L depicts the amount of cold nuclear material a product 
created in the nucleon-nucleon collision will pass through. The shading of the nuclei indicate the non-uniform density which 
should be accounted for when calculating L. 



In a nuclear environment, QVZ includes J/ip suppression via parton multiple scatterings in cold nuclear matter for 
the cc pair only, and does not consider further suppression in the hot medium created later. Multiple scattering of 
the cc pair increases the of the pair, reducingthe overlap with F{q^). The original QVZ modeled this effect with 
a constant shift linear in nuclear pathlength 0, 0| 

q^ ^q^ + Aq^ = q^ + e^L , (3) 

which was an intuitive generalization of an identical result for the nuclear broadening in photoproduction of jets given 
by twist-4 contributions 0- Surprisingly, this approach was able to reproduce not only hadron- nucleus (/i -1-^4) but 
also nucleus-nucleus {A+A) data up to SPS energies, including the "anomalous" suppression in Pb+Pb with the 
same ~ 0.2 — 0.3 GeV^/fm. The provocative result indicated negligible additional J /ip suppression from the hot 
medium in energetic heavy-ion collisions, contrary to expected signatures of a dense parton plasma. 

A more detailed analysis of coherent multiple scattering effects in the Drell-Yan process Q also found a constant 
Ag^ oc L, provided one considers at each twist only the contribution that gives the largest q^ change. The neglected 
terms, on the other hand, would affect the lower Aq^ region and therefore generate fluctuations in Aq^ at fixed L. Here 
we study fluctuations in multiple scatterings of charm quarks and antiquarks in the opposite, incoherent scattering 
limit. In that case, fluctuations in Ag^ are of the same order of magnitude as the average Ag^ for the cc pair. An 
earlier study which considered momentum transfers in three dimensions, found that fluctuations in Aq^ lead to 
much weaker suppression in A+A than in the QVZ approach. While such multiple scatterings of the heavy quarks 
are expected to have a direct impact on the J/V' Pt distributions 0,0, the QVZ and Fujii ]3l calculations do not 
provide any information about pt- In this work we focus on the interplay between the J/^/j yield and {p^) based on 
a Monte-Carlo approach. 

II. ROLE OF MOMENTUM TRANSFER FLUCTUATIONS 

In the QVZ framework, final-state quark scatterings affect the momentum distribution of cc pairs and therefore the 
J/f/' yield as 

dNj/^ = j dAq^P{Aq^;q^ - Aq^) dN,,{q^ - Aq^) = F{q^) dN,,{q^) , (4) 

where P{Aq^; q^) is the probability of accumulating a total Aq^ change in the relative pair momentum given an initial 
value q^ . A completely equivalent approach is to group scattering effects into a modified formation probability F 

dNj/^ = dA^cc(g') j dAq^P{Aq^-q^) F{q^ + Aq^) = F{q^) dN^-,{q^) , (5) 
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which has the advantage that the primary cc distribution and final state efi^ects factorize. 

The effective formation probabihty can be calculated from a model of final state interactions in the medium. We 
consider here incoherent Gaussian momentum kicks in two dimensions transverse to the beam axis, independently 
for each quark (i.e., ignore correlations between c and c), as a model of small-angle scatterings off fast partons. This 
results in momentum-space random walk probability distributions 



■K Le- 



■ exp 



(i = 1, 2 for c, c) 



(6) 



where /ct.j is the total momentum transferred, while 6^/2 is average momentum-squared transferred per unit nuclear 
pathlength. The average momentum and the average relative momentum transferred to the cc pair, kx = A;t,i + fcT,2 
and Akx = kx^i — fcT,2, are also Gaussian but with dispersions Le^. 

In the nonrelativistic limit, the accumulated change in the relative pair momentum is 



Ag^ sa Ag ^ = (g 



Ak^ 



2AkTqcos6 



(7) 



Final state scatterings can be considered at three levels of sophistication: 

i) constant shift (CS) - keep only the average of the first term, i.e., Aq^ — > Le^ 

ii) partial fluctuations (PF) - keep fiuctuations in the magnitude Afc^ but ignore the second angular term, i.e., 
Ag2 Afc|,; and 

Hi) all fluctuations (ALL) - keep fluctuations both in magnitude and direction, i.e., the full expression Q. 

In the color singlet case, the effective transition probabilities are calculable analytically in a straightforward manner: 
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where Qq = qj. — q^, q = qrf + Qg, and the last two lines are for Gaussian scattering in three dimensions. If 
fluctuations are neglected, F is strongly suppressed, exponentially with pathlength. This results in exponential J/tp 
suppression with a suppression factor Rj/^ — ex.Y>[~L e' / {2a'p)]. Fluctuations in the magnitude of Afc^ change the 
behavior to a milder power-law reduction with L, Rj/^ = 2a']p/{2a'p + L e'). The biggest effect, however, comes from 

fluctuations in the direction of Akx- These broaden F compared to the original F{q''), further weakening suppression 
effects. It is easy to understand why angle fluctuations are important. From any starting point in momentum space, 
there is always a finite probability for the random walk to get closer to the origin, thereby enhancing the formation 
probability compared to the other two approximations where any initial q' can only grow. 

Figure [2 illustrates the influence of fluctuations and confirms the above general considerations. The left plot 
shows the color singlet F(q') resulting from the different fluctuation treatments, for the original QVZ parameter set 
Nj^^ = 0.47, ap = 1.15, and with e' = 0.3 GeV^/fm and a pathlength L = 8 f m characteristic of near-central Pb + Pb 
or Au + Au collisions. The right plot shows the corresponding J/ip suppression factors (relative to the L = case) as a 
function of Le^^ for an initial cc distribution in p + p events at ^ = 200 GcV from the PYTHIA event generator 
program (the PYTHIA results can be parameterized as dNcs/dq' cx {q'f-^^^/[2m-l{f+{^q'/{l GeV') + 4,mf-^^]). 
Clearly, the largest suppression comes from the simplest, constant shift approximation (dashed-dotted line). On the 
other hand, the inclusion of fluctuations gives a much-reduced suppression (solid line). For comparison we also show 
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FIG. 2: Left plot: effective cc — > J/rp transition probability in the color-singlet channel for L = 8 x 0.3 GeV'^ (characteristic 
of near-central Pb + Pb or Au + Au collisions) with three different approaches to treat fluctuations in final-state interactions. 
Results for the constant shift (dashed-dotted) and "partial fluctuation" (dashed) approximations (see text) are compared to 
the fuU 2D Gaussian random walk result, assuming Qg = 0.5 GeV^ (solid) and the case of no final-state interactions (dotted). 
The full result is much broader in , resulting in much weaker J suppression. Right plot: suppression of the total J/ip in 
nuclear coUisions at ^/snn ~ 200 GeV as a function of Le^ for the constant shift (dashed-dotted), partial fluctuation (dashed) 
approximations, and the full 2D Gaussian random walk result (solid). The full result gives much smaller suppression than the 
two approximations. Results for the 3D random walk in 0] are also shown (dotted), for = 0.16 GeV^. 



the result for the three-dimensional Gaussian random walk considered in (dotted line), which gives a stronger 
suppression than the more realistic 2D case. The reason is that in 3D, angle fluctuations have a weaker role because 
the solid angle smOdOdcf) prefers 9 values near 7r/2, for which the angular term in Q is small. Though in the very- 
large-L limit the suppression factor behaves as 1/L and for 2D and 3D random walk respectively we found 

that in practice these limits apply poorly and give more suppression at moderate L ^ 0.5 — 5 GeV^ than even the 
constant shift approximation. All the above trends apply to the color octet case as well. 

We note that in the simple analytic calculation above we applied © with a constant Qq — 0.5 GeV^ (2D case) and 
Qq = 0.16 GeV^ (3D case), which are the average values for midrapidity J/ip-s coming from the region ^ — 3 GeV^ 
that gives the dominant contribution to the yield general, the various q components are correlated (e.g., 

Qo ^ Q^) and therefore the calculation would require dN/dq^dQ^ (2D) or dN/dq^dq^ (3D) in analytic form as an 
input. We include the q^ — Qq correlations for 2D random walk in Sec. IIIII using a Monte-Carlo approach. 

If the strength of momentum kicks per unit pathlength were an unknown parameter, much of the above discussion 
would be academic because we found that in the calculation above each approximation can reproduce the exact result 
quite well with an appropriate rescaling of e^. In particular, ejjg : e^p : e^^^ 2D '■ zd ^ ^ '■ ^-^ '■ ^-^ • ^ E3- 

However, an independent determination of from FermiLab data on dijet momentum imbalance already constrains 
£2 to - 0.2-0.5 GeVVfm 7]. 

Even if were arbitrary, it affects not only the J /ip yields but also the spectra. In particular the larger the e^, 
the higher the {p^) of the J/tp's. In the following sections we analyze the suppression- consistency in detail and 
contrast it to data from the SPS and RHIC, using a Monte-Carlo approach. For the spectra, both initial and final 
state scatterings are important. 



III. MONTE-CARLO STUDY 

In this study we utilize a Monte-Carlo Glauber model [lj| calculation, as implemented in the Heavy Ion com- 
munity which assumes a 30 mb and 42 mb nucleon-nucleon cross section for ^/snn = 17.2 and 200 GeV 
respectively, for collision geometry. We consider a nucleus-nucleus collision at a given impact parameter and collision 
energy, and for each binary collision assign a cc pair with momentum vectors from charm events in p+p collisions (at 
the same y/sj^) from PYTHIA 6.205 with the parameters used in Section^ JTJ. The pathlength depicted in Fig.^ 
L = La + Lb, is the sum of contributions by nuclcons moving towards the production point from the left and right, 
which are calculated by integrating over the nuclear density distribution and dividing by the density at the center of 
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FIG. 3: The dependence of J /ifj suppression on collision centrality for ^/snIv ~ 17.2 GeV Pb+Pb collisions for (fcy) = 0.3 GeV^. 
The filled markers indicate results from the model with changes due to re-scattering including fiuctuations. The open markers 
show results when the explicit re-scattering is replaced by an overall shift which results in the same average increase per unit 
pathlength, Aq'^/L = = 0.27 GeV^/fm. "S" (circles) and "O" (squares) refer to the color singlet and color octet transition 
probabilities (see Eqs. Q and @). 



the nucleus. E.g., for a production point (xy, zo) relative to the center of nucleus A 

La{xt,zo) = / dzp{Jz^ + x^) = ^+1^ , (9) 

where {Ncoii) is the average number of remaining binary collisions on the right and i^f^^ is the inelastic nucleon- 
nucleon cross section. Final state scatterings are modeled via a 2D Gaussian random walk as described in Section 
^except that, instead of the continuous Le^ variable, the number of scatterings is quantized in terms of the binary 
collisions, where the average momentum transfer squared per scattering is {k'^)\ i.e., in the Monte-Carlo Ncou{k^) 
is equivalent to L in the analytic calculations. We confirmed via increasing the kicks per binary collision, while 
reducing the (A;|,) by the same factor, that the discretization does not have a significant effect on our results. For 
initial state scatterings, which are meant to represent the scattering of individual incoming partons (mainly gluons) 
before the cc production, the kT is transfered to the of the cc pair with no increase in the . For simplicity, we 
assume that the (fc^) for initial state parton scatterings is the same as for the final state parton scatterings. 

From Eq. it is clear that the parameter in the Monte-Carlo is proportional to e^, the average increase in <f 
per unit nuclear pathlength. The constant of proportionality depends on the collision energy (through the inelastic 
nucleon-nucleon cross section) and also the nuclear density at the center of the colliding nuclei. For convenience, we 
tabulate in Table U several values and the corresponding e^, for Ph -f Ph collisions at ^/snn = 17.2 GeV (top 
SPS energy). 

As a consistency check, instead of our multiple scattering model, we study J/ip suppression as a function of L for 
y/SNN — 17.2 GeV Pb+Pb collisions via the (f = + e^L shift as done by QVZ. The only difference is that we use 
the PYTHIA distribution for cc pairs. The results, shown in Fig. O and labeled 'Shift', give similar suppression 
trends to those reported in I']. In particular, we observe the larger suppression as a function of L for the color octet 
transition probability F^'~'\q^) (defined in Eq. jSJ). For the Monte-Carlo studies, we use A^j/^ = 0.47 and ap = 1-2 
for F'^^Xq^) and Nj/^ = 0.485 and aF = 1 for F'^^^q^) (consistent with Sec. HHand QVZ 1]). 

Figure 13 also shows a comparison of J/ip suppression in the shifted q^ model and the fluctuation model. Just like 
for the analytic study in Section^ we keep the average q^ shift the same in both cases and therefore use {k'^) = 0.3 
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TABLE I: The = A{g^) / L values corresponding to various (fcy) parameter values for ^snn = 17.2 GeV Pb + Ph collisions. 
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TABLE IL Results for J/tp properties from the model for minimum bias collisions in y^sjvjv = 200 GeV d+Au for several values 
of input {4)- 
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GeV^ in the Monte Carlo, which corresponds to = 0.27 GeV^/fm, a shift consistent with range usedby QVZ in ^. 
We find that fluctuations give significantly less suppression. Clearly, an incoherent random walk type process cannot 
be consistent with the — + e^L shift since the most probable scattering value is zero and the low region 
therefore cannot be completely depleted. We note that our numerical results agree with the analytic ones presented 
for the singlet case in Section ITTI 

In addition, the suppression patterns for the color singlet and color octet transition probabilities, F'^^^q^) and 
^(C')(g2-j^ become nearly indistinguishable when fluctuations are included. Therefore, we only report results for the 
octet case for the remainder of this study. This is in agreement with the previous analytic calculations Q that 
incorporated momentum transfer fluctuations in three dimensions. 



IV. RESULTS 



We first use the model to study suppression and pT distributions in p{d)+Au collisions at SPS and RHIC energies. 
High statistics data are available for comparison from the NA50 experiment. In Figure 0] we show results for (fc^) 
values which reasonably reproduce either the A(p|.) = — {Pt)p+p oi' suppression as a function of the 

nuclear path-length. For p+A, a (k^) ~ 0.3 GeV^ reproduces the RpA dependence, but produces far too much A(py) 
increase. A (k^) — 0.05 GeV^, which is consistent with the A(p^) trend, under predicts the suppression for the 
largest nuclei. Even for p+A the model cannot simultaneously give satisfactory descriptions of A(p^) and RpA for 
any (fc|,) for large nuclei. We note that the relative strength of initial and final state scatterings is not tuned in 
our model, but since the initial state scatterings currently account for roughly half of the A(p|n) increase, removing 
initial state scatterings would not be sufficient to bring the model into good agreement with the data. Calculations 
for y/sNN = 200 GeV d+Au using the same {k^) values are also shown in Figure^ 

Experimental measurements of J/tp suppression RdAu and A(p^) have recently been made by the PHENIX collabo- 
ration at RHIC ITj. They report {PT)a+Au-(PT)p+p = l-77±0.37, 1.12±0.33, and -1.28±0.94 and = 1.18±0.12, 
0.79 ± 0.06, and 0.95 ± 0.1 for backward, forward, and mid rapidity respectively in minimum bias deuteron-gold reac- 
tions. Table UTI shows our Monte-Carlo results for various (fc|,) values. Although no strong statement of consistency 
can be made due to current experimental errors, the results from our model for > 0.3 GeV^ can be safely 
excluded from simultaneously satisfying both experimental RdAu and A(p^) constraints. We note that this level of 
suppression is consistent with that determined from a simple nuclear absorption Glauber model with a cross section 
^j/ip~N ~ 1 — 3 mb 0, which we have also confirmed. Future measurements of J/^ in deuteron or proton-nucleus 
reactions with sufficient statistics for centrality dependencies to be accurately determined will be important to narrow 
this parameter range. 

Despite the inability of this model to fully describe the p{d)+Au data, we feel it is useful to follow through with 
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FIG. 4: Left plots: The nuclear modification factor (Upper) and A{p|n) (Lower) for J/ip for 400 GeV/c proton beams in- 
cident on several fixed targets as a function of the mean nuclear path-length (L). The filled circles show data from NA50 
measurements |2Hl . The suppression data uses the p+p value scaled from 450 GeV/c, and the p+p reference value for {p^) 
of 1.6 (GeV/c)^ ± 10% comes from the projection of a fit to the (py) p-\-A data to L = fm. The errors for both graphs 
are dominated by the p+p reference, so they are highly correlated. Right plots: calculations for -^/sjvjv ~ 200 GeV d+Au 
collisions. Note that for the left plots each (L) comes from a different target nucleus, while for the right plots the L is for 
different production point geometries. The calculations use the F^'^\q^) transition probability. 



heavy ion results. For a high-temperature quark-gluon plasma, the incoherent final-state scattering assumption is 
more justified because the Debye length is much smaller than the mean free path ~ {gT)'^ < A ~ (g'^T)-^ 
Here we report results for Pb+Pb collisions at ^^snn = 17.2 GeV. Figure [S] shows the J/ip suppression relative to 
binary collision scaling Raa and A(p|,) as a function of centrality for various values. As (fc^) is increased, the 
fraction of surviving JM decreases and the (p^) of the surviving J/^p increases. Figure|Slalso shows data points from 
the NA50 experiment U>^. The general suppression trend is not well reproduced by any value of A value 

of (fcy) which reproduces the suppression for the most central events substantially over-predicts the increase in {p'^). 
The effect of heavy quark re-scattering in a cold nuclear medium, as described in this model, cannot explain the NA50 
data. 

Lastly, we use the model to study suppression and pT distributions in Au+Au collisions at ^/snn — 200 GeV. The 
Raa and A(p|n) results are shown for Au+Au in Fig. (HJ which again demonstrates the strong correlation between 
increased suppression, from change, and increased A(p|,) in models using this type of multiple scatterings of c and 
c quarks to account for J/ip suppression. For (fc|,) = 0.3 GeV^, roughly the largest value consistent with d+Au data, 
the fairly week suppression, Raa ~ 0.7, requires ^{Pr) ^ 5 GeV^. 

We highlight this connectedness and the inconsistency with available RHIC and SPS data for central events in 
Fig. We show that for any parameter values of (/c|,) one is constrained to be on a curve correlating suppression 
and pt- The SPS data are far from the Raa vs A(p^) curve. We note^that preliminary data from the PHENIX 
experiment also appear to be inconsistent with any points on the curve '22'] . For a given (pf^), the suppression at the 
SPS is stronger because for a lower collision energy the initial cc distribution is narrower, and therefore the yield 
is more sensitive to rescattering (i.e., shifts in q^). Although this study does not account for nuclear modifications 
of the parton distribution function (i.e., effects such as nuclear shadowing), these effects do not appear significant 
enough to change the basic conclusions of our study 0, |2^ . 
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FIG. 5: The nuclear modification factor and A{p|-) for J/?/; in y^sjvjv ~ 17.2 GeV Pb+Pb collisions as a function of the number 
of participants for several input values of {k^). The filled circles show data from NA50 measurements P.|20|. The calculations 
use the F^'^\q^) transition probability. 




FIG. 6: The nuclear modification factor and A{py) for J/tp in -^/sjvjv ~ 200 GeV Au+Au collisions as a function of the 
number of participants for several input values of (fey). The shaded area represents the range of (k^) values which are in rough 
agreement with the PHENIX d+Au data fvHi - Our calculations use the color octet F^'^\q'^) transition probability. 
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A<p^> (GeVf 

FIG. 7: The Raa vs A(j?t) curves traced by varying {k%) for events with 300-350 participants (central) for Au~\-Au and Pb+Pb 
collisions at RHIC and SPS energies respectively. The black square shows the value from NA50 measurements P. |20|I . 



V. CONCLUSIONS 

We have presented a Monte-Carlo model for calculating the effect of initial and final state re-scattering in cold 
nuclear matter on the production and suppression of J/iJj in heavy ion collisions. Although largely motivated by an 
earlier work from Qiu, Vary, and Zhang (QVZ) jj, we obtain significantly different results at a given rescattering 
"strength" (k^) due to the use of incoherent scatterings, which give considerable fluctuations in the change of relative 
cc pair momentum as opposed to the overall shift expected intuitively based on coherent re-scattering studies of jet 
photoproduction and Drell-Yan 0,0,11]. Similar observations were previously discussed by Fujii Q based on a three- 
dimensional Gaussian random walk model. Our analytic and numerical results indicate that, for a more realistic 2D 
random walk, momentum transfer fluctuations have an even stronger effect. These results highlight the importance of 
calculating the modification for a cc pair undergoing coherent multiple scattering in cold nuclear matter, including 
the first nonvanishing contribution to fluctuations in Aq^ ^ . 

In addition, we demonstrate a direct correlation between A{p^) = {Pt)a+a ~ (Pt)p+p ^^'^ ^^le J/ip suppression 
factor Raa in the incoherent multiple scattering model. We explore J/ip suppression and A(p|i) in ^/sJfN = 17.2 GeV 
Pb+Pb and ^snn = 200 GeV Au+Au collisions as a function of centrality for a wide range of {k^) parameter values, 
and find that it impossible for this class of models to simultaneously satisfy both suppression and (p^) constraints from 
available nucleus- nucleus data at SPS or RHIC. It would be interesting to compare these results with the suppression 
- (p^) correlation given by coherent multiple scatterings. 

The importance of heavy quarkonia is emp hasized by recent puzzling observations regarding open heavy flavor in 
Au+Au reactions at RHIC. Spectra [23. l25l| and azimuthal anisotropy V2{pt) = (cos 20) |2^ of "non-photonic" 
single electrons indicate a striking modification of heavy quark momentum distributions in medium. We are extending 
our Monte-Carlo calculation for J / tp by incorporating these cc pairs into the parton cascade model to thus include 
not only cold nuclear matter effects, but those of the hot medium as well. Finally, we note that our framework could 
also be applied to study nuclear effects on other heavy-flavor mesons, such as bb states, and open charm and bottom. 
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